Libraries¶
In [1]:
# libraries
import pandas as pd
import numpy as np
import datetime
import itertools
from tqdm.notebook import tqdm
from sklearn.metrics import r2_score
import lightgbm as lgb
import shap
import graphviz
# plotly stuff
import plotly.express as px
import plotly.io as pio
pio.renderers.default = "notebook"
import plotly.offline as pyo
pyo.init_notebook_mode()
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# configs
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:80% !important; }</style>"))
pd.options.display.max_rows = 50
pd.options.display.max_columns = 500
pd.options.mode.chained_assignment = None # use at your own risk..
In [110]:
pd.__version__,lgb.__version__,shap.__version__
Out[110]:
Exploratory Data Analysis¶
In [111]:
df = pd.read_csv('house_data.csv')
df.date = pd.to_datetime(df.date)
df.yr_built = pd.to_datetime(df.yr_built, format='%Y')
df.yr_renovated = pd.to_datetime(df.yr_renovated, format='%Y', errors='coerce')
df
Out[111]:
General Stats¶
In [122]:
df.info()
In [134]:
'{0:.2%}'.format(sum(df.yr_renovated.notnull()) / len(df))
Out[134]:
Missing values on yr_renovated, which is only logical. Only 4.23% completeness, precise % of expected missing values unknown at this point, will assume it's ok.
In [123]:
df.describe()
Out[123]:
Column by Column¶
In [124]:
for col in df.columns:
px.histogram(df[col]).show()
- date will be used to leave last 20% as test, to see if model is sustainable across time or suffers from decay.
- price has a lot of extreme high outliers. We might be better of capping the training set up to X value and add to model deployment a caveat that states that anything at X is intented to be read as X or above.
- In theory, performance metrics would be almost the same, but in practice this preprocessing could help for better model implementation.
- bedrooms/bathrooms have the odd 0 or extreme value outliers, which can be handled just fine with the tree based approach LightGBM.
- sqft_living/lot/above/basement, same logic applies.
- zipcode can be used as region identifier. Must be category-encoded for LightGBM first.
- lat and long are great as they are to serve as geographical mapping for LightGBM.
- sqft_living15/sqft_lot15 no info
- Unmentioned colums look fine.
sqft_15¶
In [97]:
px.scatter(df.sample(1000), x='sqft_living15', y='sqft_living', trendline='ols')
In [96]:
px.scatter(df.sample(1000), x='sqft_lot15', y='sqft_lot', trendline='ols')
- There's a clear trend between the two, but I have no information at this point as to what 15 means.
- Could be a good idea to eventually code brute force search to weed out underperfroming features.
Anomalies Checks¶
renovated <= built ?¶
In [90]:
df[df.yr_renovated <= df.yr_built].sort_values('yr_renovated')
Out[90]:
Looks fine, no houses renovated before they were built.
sqft_living != sqft_above + sqft_basement ?¶
In [101]:
df[df.sqft_living != df.sqft_above + df.sqft_basement]
Out[101]:
OK, no anomaly on here.
sqft_living > sqft_lot ?¶
In [103]:
df[df.sqft_living > df.sqft_lot]
Out[103]:
In [104]:
px.scatter(df[df.sqft_living > df.sqft_lot], x='sqft_living', y='sqft_lot')
Since living = above + basement checks out, I will assume these datapoints are wrong and clean them as null.
Data Preparation¶
Analysis¶
Nullify Values¶
- df[df.sqft_living > df.sqft_lot] as explained above.
Feature Engineering¶
If we simply use the year for built/renovated columns, this will be progressively outdated in the future. We must use the time between that data and the data capture date. If the divisor is zero or NaN and the calculation returns NaN, LightGBM handles just fine.
- yr_built_ago = years(date - yr_built)
- yr_renovated_ago = years(date - yr_renovated)
- yr_ren_blt_ago_dif = yr_renovated_ago - yr_built_ago
- yr_ren_blt_ago_rat = yr_renovated_ago / yr_built_ago
Since LightGBM is a tree based approach which suffers from both additive and linear limitations, differences and ratios to represent seemingly important relationships can help.
- bath_bed_sum = bathrooms + bedrooms # sum of rooms
- bath_bed_rat = bathrooms + bedrooms # bathrooms per bedroom
- bed_floor_rat = bedrooms / floors # rooms per floor
- bath_bed_sum_sqft_living_rat = sqft_living / bath_bed_sum # sqft per room
- floors_sqft_living_rat = sqft_living / floors # sqft per floor
- sqft_basement_rat = sqft_basement / sqft_living # relative size of basement to living area
- sqft_living_rat = sqft_living / sqft_lot # relative size of living to lot area
- sqft_living_15_rat = sqft_living / sqft_living_15 # we're comparing something but don't know what, hopefully this helps, should check
- sqft_lot_15_rat = sqft_lot / sqft_lot_15 # we're comparing something but don't know what, hopefully this helps, should check
Scores generated from relevant, to overcome additive limitation of LightGBM which can't straightforwardly combine both.
- water_view_score = waterfront * 5 + view # waterfront gets 5x weight to equally weight to view
- condition_grade_score = condition * 2 + grade # condition gets 2x weight which roughly puts it on the same scale as grade
Inputs Exclusion¶
- id, looks like a random id number with no value.
- yr_built/renovated, they are used as timedelta in F.E.
Functions¶
In [2]:
def nullify(df):
df.sqft_lot = np.where(df.sqft_living > df.sqft_lot, np.nan, df.sqft_lot)
return df
def feng(df):
# categorical conversion
df.zipcode = pd.Categorical(df.zipcode)
# built and renovated
df['yr_built_ago'] = (df.date - df.yr_built).dt.days / 365.25
df['yr_renovated_ago'] = (df.date - df.yr_renovated).dt.days / 365.25
df['yr_ren_blt_ago_dif'] = df.yr_renovated_ago - df.yr_built_ago
df['yr_ren_blt_ago_rat'] = df.yr_renovated_ago / df.yr_built_ago
# rooms
df['bath_bed_sum'] = df.bathrooms + df.bedrooms
df['bath_bed_rat'] = df.bathrooms / df.bedrooms
df['bed_floor_rat'] = df.bedrooms / df.floors
df['bath_bed_sum_sqft_living_rat'] = df.sqft_living / df.bath_bed_sum
df['floors_sqft_living_rat'] = df.sqft_living / df.floors
# sqft
df['sqft_basement_rat'] = df.sqft_basement / df.sqft_living
df['sqft_living_rat'] = df.sqft_living / df.sqft_lot
df['sqft_living15_rat'] = df.sqft_living / df.sqft_living15
df['sqft_lot15_rat'] = df.sqft_lot / df.sqft_lot15
# extra points
df['water_view_score'] = df.waterfront * 5 + df.view
df['condition_grade_score'] = df.condition * 2 + df.grade
return df
def exclusions(df):
del df['id']
del df['yr_built']
del df['yr_renovated']
return df
In [26]:
def r2_metric(preds, train_data):
labels = train_data.get_label()
r2 = r2_score(labels,preds)
eval_name = 'r2'
eval_result = r2
is_higher_better = True
return eval_name,eval_result,is_higher_better
def lgbm_model_xv(df, inputs, output, params_grid):
params_df = pd.DataFrame()
tr_ds = lgb.Dataset(df[inputs], label=df[output])
for params in tqdm(params_grid):
cv = lgb.cv(
params = params,
train_set = tr_ds,
num_boost_round = 10000,
early_stopping_rounds = 25,
verbose_eval = False,
stratified = False,
eval_train_metric = True,
feval = r2_metric,
seed = 7,
)
cv = pd.DataFrame(cv).tail(1)
params = pd.DataFrame(params, index=[0])
params['num_boost_round'] = cv.index.values[0] + 1
params['perfo'] = cv.iat[0,2]
params['perfo_train'] = cv.iat[0,0]
params['perfo_vw'] = 2 * params['perfo'] - params['perfo_train'] # this is variance-weighted performance
params = pd.concat([params,cv.reset_index(drop=True)], axis=1)
params_df = params_df.append(params)
return params_df
def get_best_params(params_df, ascending):
params_df = params_df.sort_values('perfo_vw', ascending=ascending).tail(1)
perfo = params_df['perfo'][0]
perfo_vw = params_df['perfo_vw'][0]
params_df = params_df.iloc[:,:params_df.columns.get_loc('perfo')] # get all columns until perfo data comes in
params = params_df.to_dict('index')[0]
return params,perfo,perfo_vw
def train_lgb(df, inputs, params):
rounds = params[0]['num_boost_round']
params[0].pop('num_boost_round', None)
tr_ds = lgb.Dataset(df[inputs], label=df[output])
booster = lgb.train(
params = params[0],
train_set = tr_ds,
num_boost_round = rounds
)
return booster
In [112]:
# read input
df = pd.read_csv('house_data.csv')
df.date = pd.to_datetime(df.date)
df.yr_built = pd.to_datetime(df.yr_built, format='%Y')
df.yr_renovated = pd.to_datetime(df.yr_renovated, format='%Y', errors='coerce')
# general data preparation
df = nullify(df)
df = feng(df)
df = exclusions(df)
# testing flag
df = df.sort_values('date')
df['testing'] = 1
df.testing = df.testing.cumsum()
df.testing = np.where(df.testing < df.testing.max() * 0.8,0,1)
# define inputs and output column
inputs = list(set(df.columns) - set(['date','price','testing']))
output = 'price'
# parameters to explore
params_grid = {
'metric': ['r2_metric'],
'objective': ['mse','mae'], # mae penalizes outliers less, let's see who wins
'num_leaves': [4,8,16,32],
'learning_rate': [0.1],
'feature_fraction': [0.5,1],
'bagging_fraction': [0.5,1],
'bagging_freq': [1],
'verbosity':[-1],
'seed':[7],
}
params_grid = [dict(zip(params_grid.keys(), v)) for v in itertools.product(*params_grid.values())]
Model Training¶
In [18]:
# run parameter search
params_df = lgbm_model_xv(df[df.testing==0],inputs,output,params_grid)
In [25]:
params_df = params_df.sort_values('perfo_vw', ascending=False).reset_index(drop=True)
params_df
Out[25]:
- The best model according to perfo_vw is shown first.
- This metric, perfo_vw, is a measure of performance I created to value less-variance models over simply selecting lowest bias.
- I'll select model index 4 here, sacrificing a little bit of variance in order to ramp up the bias almost to 0.9.
Let's now use these parameters to train the final model to be tested on the test, out of time sample.
In [56]:
# train chosen model
params = pd.DataFrame(params_df.iloc[4]).transpose().reset_index(drop=True)
params = get_best_params(params, True)
booster = train_lgb(df[df.testing==0], inputs, params)
Model Testing¶
In [96]:
df_test = df[df.testing==1].copy().reset_index(drop=True)
df_test['price_yhat'] = booster.predict(df_test[inputs])
In [113]:
perfo_tr = params[1]
perfo_te = r2_score(df_test.price,df_test.price_yhat)
print('Training R2:',round(perfo_tr,4))
print('Testing R2:',round(perfo_te,4))
print('Loss:......:','{0:.2%}'.format(1 - perfo_te / perfo_tr))
2.15% loss for out of time testing looks acceptable for me, although more exploration could be done.
In [99]:
mape = df_test[['price','price_yhat']]
mape['percentage_error'] = abs(mape.price - mape.price_yhat) / mape.price
'{0:.2%}'.format(mape.percentage_error.mean())
Out[99]:
12.77% MAPE also looks acceptable, although further investigation could be in regards to in what price range is the most error coming from.
In [100]:
px.scatter(df_test, x='price', y='price_yhat', trendline='ols')
Residuals visual analysis shows that the model is doing pretty well overall, as seen before with R2 and MAPE.
Model Discovery¶
Let's see how we can explain our model. Since it's a tree-based algorithm, we can start be taking a look at the first tree.
In [101]:
lgb.create_tree_digraph(booster, tree_index=0)
Out[101]:
A simple tree-based algorithm simply takes the features given and tries to split them to minimize the chosen loss function and in turn maximize the chosen performance metric. This tree looks really nice and with some effort we could extract information, but our algorithm is quite complex and actually has 791 trees which are aggregated to compute the final prediction. We can't go through all of them. Therefore, we will use SHAP, a package with a game-theory approach to explaining complex, black box models such as this one.
In [116]:
# shap init
shap.initjs()
explainer = shap.TreeExplainer(booster)
shap_values = explainer.shap_values(df_test[inputs])
We can see at a record level which features are driving changes up or down in price, which has a base value of its own mean at 532k.
In [128]:
shap_plot = shap.force_plot(explainer.expected_value, shap_values[0], df_test[inputs].iloc[0])
shap.save_html('shap_force_single.html', shap_plot)
shap_plot
Out[128]:
- West side longitute lowers value. Indeed, in general, left side is cheaper.
- But, north side latitude drives value up, as was also seen on the map. You can review this below on the scatterplot.
- A particular zipcode drives value way up.
Overall, all what we are seeing makes sense to what we expect, but this approach allows us to see exactly what are the score components and how much they affect the final score. In this way, a user analyst from the real state can take a look at this 'force-plot' and make the final decision himself, by understanding why the model is settling at a particular value and what information it could be missing.
In [118]:
px.scatter(df_test, x='long', y='lat', color='price_yhat')
What about overall feature impact on price? Let's take a look at a much more complex visualization:
In [129]:
shap_plot = shap.force_plot(explainer.expected_value, shap_values[:1000], df_test[inputs].iloc[:1000,:])
shap.save_html('shap_force_multiple.html', shap_plot)
shap_plot
Out[129]:
You can see here lots of things, such as:
- lat effects by lat: there's not a linear relationship. As seen on the map, lat up to mid-north drives price down, then up, but towards the north it again goes down. This is a classic non-linear relationship that a tree algorithm can find.
- lat effect by sample order by output value: as this is a complex multivariate model, some high-price houses can have drive-down lat and viceversa. These cases surely compensate with other opposite effects to get to final value.
Finally, let's take a look at overall feature value based on impact on model output.
In [131]:
shap.summary_plot(shap_values, df_test[inputs], max_display=len(inputs))
Our condition+grade and water+view scores made it pretty far up, good job!
Proposed Improvements¶
Feature Extraction / Engineering¶
- Since zipcode and location is so important, we could get other geodata such as neighborhood.
- Since grade is so important, we could see if the data that generates such aggregated grade is available and use that higher detail information instead.
- Since the sqft features are so important, we could try to get more data in relation to that, such as sqft of the rooms.
- Actually I also missed the yard sqft calculation, which is a differentation which would help LightGBM because of it's additive limitation.
- Shores distances.
- Get arrays of geopoints to delimit each of the relevant lakes and water areas, such as Lake Washington and the estuary.
- For each delimited area, compute minimum distance between each house's geopoints and the area's array of geopoints.
- Add to model and see if there's a performance increase.
- Lot height.
- https://ngmdb.usgs.gov/ngmdb/ngmdb_home.html from here topographical maps can be obtained and used to calculate lot height.
- We should have a discovery session with the real state experts and see what features we are missing to obtain near-perfect performance.
Machine Learning¶
- Search for a wider parameter space on LightGBM and lower learning rate.
- Try more algorithms and ensemble algorithms of multiple learners.
- Try bruce-force feature subspace search. Sometimes features can be selected by the tree algorithm but actually are detrimental to XV performance.